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Boundary driven diffusive systems describe a broad range of transport phenomena. We study 
large deviations of the density profile in these systems, using numerical and analytical methods. 
We find that the large deviation may be non-difFerentiable, a phenomenon that is unique to non- 
equilibrium systems, and discuss the types of models which display such singularities. The structure 
of these singularities is found to generically be a cusp, which can be described by a Landau free 
energy or, equivalently, by catastrophe theory. Connections with analogous results in systems with 
finite-dimensional phase spaces are drawn. 
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I. INTRODUCTION AND FRAMEWORK 

The dynamics in many systems of physical interest are 
described by a field p{x,t), with large-scale conserving 
diffusive behavior and noise. For example, p{x, t) could 
describe the density of diffusing particles, the local tem- 
perature in a heat transport experiment, or any other 
field which behaves diffusively. For such systems, when 
the interactions are short range, it is accepted that 
the large-scale behavior of the current obeys Fick's- (or 
Ohm's- or Fourier's-) law with noise. Here our interest is 
in transport experiments, where the system is attached 
to reservoirs, whose effect is to fix the value of p at the 
boundaries, resulting in a net current flowing down the 
gradient. 

For such systems the density p{x,t) and the current 
J (x, t) satisfy the conservation relation 



dtp + d^J^O 



(1) 



where 



J=-D{p {x, t)) d^p [x, t) + v/fT (p (x, t))77 {x, t) . (2) 

Here D{p{x,t)) is a density-dependent diffusivity func- 
tion, and (j{p{x,t)) controls the amplitude of the 
white noise •q{x,t), which satisfies {rj{x,t)) — and 
(77 {x, t) rj {x', t')) = {x - x') 6{t- t'). At temper- 

atures well-above any phase-transition, which we study 
here, D (p) and a (p) are smooth functions, and D > 0. 
For simplicity we consider one dimension, where the dis- 
tance is rescaled by the system size N, so that < 
a; < 1, and time is rescaled by N'^. The small 
term in the noise is a direct consequence of this coarse- 
graining. D [p) and a (p) are related via a fluctuation- 
dissipation relation, which for particle systems reads 
a (p) = 2kBTp^n (p) D (p) where k (p) is the compress- 
ibility f]\. The system is attached to reservoirs at the 
boundaries x = 0, 1, which act as boundary conditions 
(BCs), p{x = 0,t) = PL and p{x = l,t) = pn- If Pl ^ PR 
a current is induced through the system, driving it out 
of equilibrium. For applications of Eq. ([2]) to transport 
phenomena, including electronic systems, ionic conduc- 
tors, and heat conduction, see for example [1,[1,0|- 



It is natural to ask for the probability of a density pro- 
file pf (x) in the steady-state. It is known that the prob- 
ability distribution assumes the form P[pf] ~ e^^'^[''^l, 
where 4>\pf] is known as the large deviation functional 
(LDF), and the N prefactor is due to the small noise. 
P[pf] is the subject of this work. As seen, out of equi- 
librium the LDF [p/] Jjlays the role of the free-energy 
density in equilibrium 1]. It is by now well established 
that, in contrast to equilibrium, out of equilibrium long 
range correlations build-up [1, Q and the LDF is non- 
local Moreover, there is now a general framework 
for calculating (/>[p/]. As detailed below it involves find- 
ing the most probable history leading to pf. In spite of 
this, the properties of 4>[Pf] remain poorly understood. 
Much of what is known is based on a handful of exact 
solutions for specific models p^| - [T2| , and numerical tech- 
niques [isj . In particular, it was recently realized that 
can be non-differentiable pll ]. 

Here we discuss in detail the occurrence and structure 
of such singular behavior in the class of models defined 
above. We refer to it as a Large Deviation Singularity 
(LDS). This is very different from the equilibrium case, 
where smooth dynamics (i.e. when D (p) and a (p) are 
smooth and D > 0) lead to a smooth LDF [p/]. 

The general occurrence of non-differentiable LDFs in 
low-noise Langevin equations was discovered by Graham 
and Tel LDSs were consequently widely discussed 
in the literature [l5l - [l9| . demonstrated in experiments 
[Tsl [T6j . and shown to affect quantities such as barrier 
crossing rates [13] ■ In addition to these, and more closely 
related to the present work, an LDS was proven to exist 
in the asymmetric exclusion process [20| , a specific model 
of diffusing particles, where (unlike in the present paper) 
the particles are also subject to an driving field in the 
bulk. This is perhaps the first known microscopic lattice- 
gas model for which the continuum limit was proven to 
feature an LDS. The proof hinges on the exact solvabil- 
ity of the model, and it is not clear what other models 
of that family will show this behavior. The general con- 
ditions for LDSs to occur in fields, and the structure of 
the singularities remains largely unknown. 

In this paper, we achieve the following. 

(1) Existence of non-differentiable LDFs for boundary- 
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driven systems - we show that in some boundary-driven 
diffusive models the LDF is non-differentiable. This in- 
cludes an exactly solvable model, where D = 1 and 
a = + 1, and a boundary-driven Ising model, with 
conserving dynamics in the bulk. The phenomenon is 
general and robust, and expected to be found in models 
where a {p) / D (p) has a minimum which is deep enough. 
The profiles pf (x) at which the derivative S(f)/6p is dis- 
continuous are found to have a typical shape, as shown 
in Fig. [2Jb) and SI The jump in the derivative is due to 
a change in the form of the most probable history p (x, t) 
leading up to pf (x). It stems from the existence of re- 
gions in the space of /?/ (x) where multiple locally mini- 
mizing histories lead to a single pf (x). A short account 
of these results was given in j2ll |. 

(2) Structure of singularities in phase-space - we study 
the singular structures in phase-space. Two-dimensional 
cross-sections, as shown in Fig. ^c) and [SJa), are il- 
luminating. They show the regions in the pf (x) space 
with multiple locally minimizing histories. The bound- 
aries of this region are known as caustics. The points 
where 5(j)/5p jumps occur when the globally minimiz- 
ing history changes. These points form the transition 
line. The transition line and caustics end at a sin- 
gle point, analogous to a second-order phase transition. 
We show that the structure near this point is similar 
to the one described by a Landau mean-field theory, or 
by a cusp singularity in catastrophe theory. One out- 
come of this theory is the prediction that at the second- 
order-like point the probability P[pf] scales in a non- 
analytic way with the system size N . Specifically, instead 
of the expected series lnP[p/] = —N(t>[pf] -\- 0{N^^, 
the series will have an additional logarithmic correction 
In P[p/] = ~N(t)[pf] -1- l/41niV -I- O {N°). The prefactor 
1 /4 is universal, depending only on the symmetries of the 
systems. 

(3) Relation to finite dimensional theory - we show 
that much of the physics can be understood by intro- 
ducing simple toy models with as low as two degrees of 
freedom. 

We stress that the singularities discussed here are dif- 
ferent in nature from those found for global quantities 
such as the current [22| - [25t . In those case the probability 
of a configuration can be smooth in phase-space, but the 
optimizing configuration can change abruptly. 

The paper is arranged as follows: In Sec. |n]we present 
an example of an LDS in a model with a single degree of 
freedom, as well as the background to the general theory. 
In Sec. Illll we demonstrate the existence of LDSs in mod- 
els of the family discussed here. We introduce the two 
example models which are studied throughout the paper. 
We show, analytically for one model and numerically for 
the other, that LDSs do indeed exist, and indicate where 
and under what conditions they are expected. In Sec. 
ITVlwe study the structure of the region, and the effect of 
this structure on the dependence of the probability P[pf] 
on the system size. In Sec. |V]we introduce a model with 
just two phase-space dimensions which, as we show, cap- 
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FIG. 1. Simple model with singular LDF. The gray curve is 
U {x) — —fox — V (x), and black curve is the LDF (x). The 
dashed arrows show the most probable path for a particle from 
A, the local minimum of the potential, to a point between B 
and C. 0' (x) is discontinuous at point C. 



tures much of the behavior of the full infinite-dimensional 
field model. In Sec. IVII we conclude and discuss future 
research directions. 



II. LDSS IN SIMPLE MODELS 

Before discussing LDSs in the model defined above, 
we recall the simplest example of such a phenomenon, 
which occurs for a single particle movin g on a ring. As 
originally discussed by Graham and Tel when such 
a system is driven out of equilibrium the gradient of the 
LDF becomes discontinuous. It is instructive to see how 
this singularity arises, despite the fact that many of the 
features are different from the singularities discussed in 
this paper. 

Consider a particle moving on a one dimensional ring 
X G [0, 1], subject to the Langevin equation 

rj nf 

- = fo-V'{x) + V~ev{t) , 

where (77 (t) rj {t')) = 26 {t — t'), /q is a constant, V (x) a 
periodic function on the ring, and e is a small number 
which plays an analogous role to in Eq. ([2]). For 
/o 7^ 0, the total force F (x) — fo~V' (x) is not derivable 
from a potential. It is useful to define the integral U (x) — 
F {x') dx' for X e [0, 1] which is no longer periodic in 
X. Consider the case with U (x) shown in Fig. [TJ 

As the noise is small (because of the e prefactor), the 
system spends most of its time near point A. We now 
want to evaluate the probability P (x) ~ "^^^^ of 

reaching any other point, in the low noise limit to leading 
order in e~^. The point B is represented both by 2: = 
and X = 1. However, due to the bias force it is easier to 
reach B by moving to the right. Therefore the probability 
of reaching point B is given by the Arrenius factor 

p _ ^-e-'[U{l)-U{A)] ^ g-e-i0(i3) ^ 

To see why cf) (x) is singular, note that once the particle 
has reached this point it may "roll-down" to reach all 
points. Therefore the probability of being between B and 
C is equal to this order. Only below C is it preferable 
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to move from point A to the left, and the probabihty 
changes again. Thus for /q ^ one obtains a plateau 
and a discontinuity in 0' (x) at point C. 

While the plateau is a rather specific feature of this and 
similar examples, the existence of a discontinuity in the 
derivative 0' (x) is a common feature of non-equilibrium 
low-noise systems. It results from a competition between 
different trajectories which lead to the same final point. 
Here, these are trajectories moving to the right and left. 

In the context of boundary-driven diffusive systems, 
the LDFs of all models which were previously studied 
exhibited a smooth LDF. This raises the question of 
whether, and for what models of this family, will LDSs ex- 
ist. In addition, it is interesting to understand the struc- 
ture of these singularities, and whether they are similar 
to what is known for models in a finite-dimensional space, 
where the structure can be understood using mean-field, 
or catastrophe theory. 



A. Background theory 

We now outline the theoretical tools used below. The 
average, or most probable density profile for the system, 
p, is obtained by solving dx [D {p) dxp] — 0, with p (0) — 
PL and p(l) = PR at the boundaries. As the noise is 
small, the system spends most of its time close to p. In 
order to find the probability of any profile pf (x), one 
must therefore compute the probability of reaching pf, 
starting from p in the distant past. 

The probability density of a history {p (a;, t) , J (x, i)} 



during time —oo < t < is P - 
S is given by [S-SnnTEil 



where the action 



5= r dt [\^ [JM+D{pM)dxP{x,t)f 
J-oo 2aipix,t)) 

. . (3) 

The probability density P [pf] of reaching p/ is then given 
by the path integral 



P[pf]^ J Dp J 



DJ^-NSlp,J] 



taken over histories satisfying dtp + dxJ — 0, with initial 
and final conditions p{x,t ^ —oo) = p (x), p (x, t = 0) = 
Pf{x), and boundary conditions p(x = 0,t) = p^ and 
p {x — l^t) = pif. For large TV a saddle-point approxima- 
tion gives P ^ e"^'^^''^] with [pf] = infp.j S, where the 
infimum is over all allowed histories. 

In equilibrium (i.e. when p^ — Pr), the steady-state 
probability of a density profile pf (x) is easy to obtain - 
the LDF (j) [p f] is then given by the free-energy which is 
local in Pf, 0[/O/] = J f {p f {x) , p) dx, where 



dpi 



dp2 — ^-Y- 



(4) 



Note that in this case p is constant, p = Pl = Pr- By con- 
trast, the steady-state probability distribution away from 



equilibrium is notoriously hard to compute, and very dif- 
ferent from the naive guess (f>[pf] — J f {p f {x) , p (x)) dx, 
now with space dependent p (x). In fact, as stated above, 
(j)[pf] is non-local. 



III. EXISTENCE OF LDS 

As is clear from Eq. LDSs cannot exist in equi- 
librium if D (p) and a (p) are smooth, positive functions. 
We now turn to discuss non-equilibrium cases where they 
can exist. In previously studied exactly solvable non- 
equilibrium models [13, El] , the action S in Eq. (jS]) has 
a single local minimal history leading to pf, and (f) [pf] is 
then a smooth functional. However, this need not always 
be the case, and there can be multiple local minima to 
S to the same pf. For example, in the model of a single 
particle described in Sec. [Hi these are trajectories corre- 
sponding to the particle moving left or right on the ring. 
When this occurs, the global minimum can switch be- 
tween the different local minima (as it does at point C in 
Fig. [T]). For fields, this is accompanied by a jump in the 
functional derivative of the large-deviation Scp/Spf. This 
is reminiscent of the mechanism for a first-order phase 
transition in equilibrium. 

It is unclear which models display an LDS. It is known 
that models which feature more than one fixed point of 
the zero- noise dynamics generically display LDSs p^ . 
In the models we study here, however, p is the only 
zero-noise fixed-point. Therefore, the absence of LDSs 
in previously studies models of this class is not surpris- 
ing^,'l^. 

In this section we discuss the existence of LDSs in two 
models. One of the models originates by taking the con- 
tinuum limit of a "microscopic" lattice gas model, the 
driven Ising model. The other has the advantage of be- 
ing exactly solvable. As explained in Sec. Ill Al a model 
is defined by two functions, the diffusivity D (p) and the 
noise strength a {p). These are shown in Fig. [H^a) for a 
specific set of parameters of the driven Ising model, and 
Hfb) for the analytically solvable model, referred to as 
the quadratic-cr (QS) model. We first define the models, 
and then discuss their properties. 



A. Definition of models 

Here we define two models for which we demonstrate 
the existence of LDSs. As a common feature, both mod- 
els display a pronounced dip in the function a {p) / D (p) . 

1. The boundary- driven Ising (BDI) model 

This is a lattice gas with on-site exclusion and nearest- 
neighbor interaction. It corresponds to the Katz- 
Lebowitz-Spohn model [13] with zero bulk bias. The 
model is defined on a Id lattice with sites i = l,..,iV, 
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FIG. 2. Model definitions. The functions g (p) and 73 (p) for 
(a) the BDI model, and (b) the QS model. 



each of which can be either occupied ("1") or empty 
( "0" ) . The model depends on two rate parameters 5 and 
e. The jump rate from site i to site i + 1 depends on the 
occupation at sites i — 1 to z + 2 according to the following 
rules: 

0100 0010, 1101 1011 , 
1100 1010, 1010 1100 , 



form will be used throughout the text. Note that the 
BCs of the density map accordingly. 

The QS model has the advantage that it is analytically 
tractable [l^: the LDF is given by 4>[pf] — miiKpext, 
where (j^ext are extremal values of the action given by 



dx <f{pf {x),g{x)) - In 



P'{x) 



(5) 



Here / (p, g) is defined in Eq. Q and g {x) is an auxiliary 
function satisfying the differential equation 







g [x) - pf (x) g" {x) 



(6) 



with BCs g (0) — pL, and 5(1) = pr. Note that a.s D ~ 
1, the most probable configuration p{x) is linear, with 
p(0) = PL and p(l) = PR. Each of the solutions of Eq. 
(Et , when used in Eq. ^ , gives (pext of an extremal path 



B. The use of cross-sections 



and their spatially inverted counterparts with identical 
rates. 

For equilibrium BCs, e.g., periodic BCs, the dynamics 
admits an Ising probability distribution P oc cxp {—fiE) 
with 



S = ^ (1 - 2n,) (1 - 2n,+i) + M^I 



This energy describes nearest neighbor interactions, and 
a chemical potential term. /3 is related to e by exp (4/3) = 
(1 + e) / (1 — e), and /x fixes the average density. The pa- 
rameter 5 does not affect the stationary state, but does 
enter into the dynamical behavior of the model. For each 
parameter set (e, 5) one can write implicit analytic equa- 
tions for D (p) , (J (p) which can then be inverted numer- 
ically. The calculation is described in Appendix [A] Fig. 
IDJa) shows D (p) and a [p) for (e, 5) = (0.05, 0.995). As 
can be seen, D (p) is peaked and a (p) has a local mini- 
mum around p = 1/2. This will be a key feature of the 
model. 



The quadratric-(T ( QS) model 



The model is defined by constant D and a (p) — 
c{p— b)^ + a, with a, c > 0, so that a (p) is a parabola 
clear above the axis. Upon the rescaling 

/a fa 
-p + b , J DJ-J , 

x-^ X , t/D , S cS/D 

the model can be brought to a standard form defined by 
D = 1 and cr (p) = p2 -h 1, see Fig. [^b). This standard 



Below we demonstrate the existence of LDSs in the 
models defined above. As the phase-space is infinite di- 
mensional, the structure of is hard to visualize. For 
many purposes it is sufficient to consider two-dimensional 
cross-sections of the infinite-dimensional phase-space. 

To this end, in most of what follows, out of the phase- 
space of final profiles p / {x) we restrict ourselves to those 
parametrized by just two variables, of the form 



p / {x) = p [x] + ai sin t:x + a2 sin 2-kx 



(7) 



This is a cross-section in the phase-space of final states. 
(Note that in Appendix [B] we use a different form.) It 
will be more convenient to parametrize these profiles us- 
ing p/ (1/3) and p/(2/3) instead of Q!i,a2. We stress 
that this choice is rather arbitrary and that the singular- 
ity described occupies a space of co-dimension 1 in the 
infinite-dimensional phase space. 

In order to visualize trajectories p {x, t) leading to 
Pf (x), we plot p{x = 2/3, i) against p{x = 1/3, t). Note 
that here we do not constrain p (x, t) at intermediate 
times to be of the form in Eq. (O . 



C. Non-unique path minimizers and the LDS 

As we now show, in both the BDI and the QS models, 
there are certain states p/ for which there exists more 
than a single history p (x, t) that extremalizes the action 
in Eq. ([3]). In order to find multiple extremal solutions 
we use different techniques, depending on the model. 

In the QS model we look for solutions to the differential 
equation ([5]). These are found using a shooting method 
[29|, in which Eq. ([5]) is integrated from a; = to a; = 1, 
with initial conditions i? (0) = pi, and 5' (0) — c. The 
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FIG. 3. QS cusp, (a) A profile p/ (solid line) for which Eq. 
((5]) has a single solution (dashed line), (b) A profile p/ for 
which Eq. ^ has three solutions. Pf of panel (a) is shown 
for comparison (dotted line), (c) A cross-section - the density 
pf {x = 2/3) vs. pf {x = 2/3), for configurations of the form 
given by Eq. 0. The region of multiple solutions (gray), and 
the switching line (dashed line). Crosses denote the locations 
in this cross section of the profiles shown in panels (a,b). 



values of c are scanned systematically to find all solutions 
where g {!) — pR. In this way all solutions of Eq. ([5]) are 
obtained. 

Using final profiles p/ of the form of Eq. ([7]) with 
PL — —3, PR = 3, we find two distinct behaviors. For 
final profiles which lie in the white region of Fig. [3jc) 
we obtain a single solution to Eq. ([6]), as illustrated in 
Fig. Ela). In contrast, for final profiles in the gray region 
Fig. ^c) we find three solutions, see Fig. EJb). Of these 
three solutions, two correspond to local minima of the 
action Eq. ^ and one to a saddle-point. Among the 
two minima, one is lower than the other except along 
the switching line, marked by a dashed line in Fig. |3^c), 
where they are equal. On this line the global minimum 
switches from one local minimum to the other. This leads 
to a jump in the gradient of the LDF S(f)/6pf across the 
line. 

The phase diagram, shown in Fig. [3l[c), is reminis- 
cent of that obtained from a Landau free-energy. In this 
analogy, the gray region corresponds to the free-energy 
having two local minima, one metastable. The bound- 
aries of the gray region are then the spinodal lines (where 
the metastable minimum disappears), and the switch- 
ing line corresponds to a first-order transition (where 
the metastable and stable minima exchange roles). The 
switching line terminates at a point analogous to a crit- 
ical point. We examine this issue in detail below, and 
show that a universal behavior emerges. 

It is natural to ask which BCs admit profiles Pf{x) 
with multiple minimizing solutions. In the case of the QS 
model, we can in fact show that: For any BCs pL ^ Pb., 
there exists a profile pf (x) for which Eq. (0) is satis- 
fied by more than one solution. The proof is given in 
Appendix |B] This is interesting since it implies that in 
this model even the smallest deviation of the BCs from 
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FIG. 4. LDS in the BDI model. Two locally minimizing 
histories leading to the same p/, plotted at different times. 




FIG. 5. LDS in the BDI model, (a) The time evolution of 
p(l/3, t) vs. p(2/3,t) is plotted for the histories of Fig. |4l 
(b) The switching line (dashed line), together with lines of 
equal (solid lines). 



equilibrium leads to the existence of LDSs. The closer 
the BCs are to equilibrium, the further the states pf are 
from p before multiple solutions exist. 

We now turn to the BDI model. Here no analytical 
solution is known, and we solve for local minimizers of 
the action S, using the numerical technique described in 
[H, mi- Again, for 7^ pR we find final configurations 
Pf with multiple minimizing solutions. Fig. Ogives an 
example of such a pj. The two paths leading to this con- 
figuration are also shown in Fig. [Sfa), where we plot the 
values of p (x = 2/3, t) against p{x — 1/3, t) of the same 
histories. In Fig. El^a) we also plot the numerically ob- 
tained region in the pf (1/3) and p/ (2/3) plane for which 
multiple histories are found, as well as the switching line. 
The jump in the gradient (50/(5py is clear in Fig. [2Ub), 
which depicts lines of equal 0. 

The LDSs in the two models have many features in 
common. The phase diagrams in Fig. [3fc) and[5t^a) have 
a similar structure, with the singularities appearing for 
similar final profiles p/. There is one important differ- 
ence: In contrast to the QS model, in the BDI model 
a finite difference of the boundary conditions pR — p^ is 
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needed in order for an LDS to exist. In both models, gen- 
erally we find (data not shown) that as the value pn — p^ 
is decreased, the region with multiple solutions is pushed 
away from p. However, in contrast to the QS model 
where p is unbounded, in the BDI model p is bounded 
(0 < /O < 1). Hence below some threshold value, no LDS 
is found in the BDI model. Similarly, by tuning e and 6 
in the BDI model, D and a can be continuously varied 
from the simple symmetric exclusion model with D ~ 1 
and a = 2p(l — p), for which the LDF is known to 
be smooth, to the model discussed above. The singular- 
ity appears when the dip in cr (p) / D (p) is deep enough 
(data not shown). 

To summarize, in both models we find LDSs when the 
function cr (p) / D (p) has a (deep enough) local minimum. 
Numerical experiments indicate that this is indeed, more 
generally, the requirement. Recall that by fluctuation- 
dissipation, the ratio is related to the compressibility 
a [p) /D (p) = 2kBTp'^K (p). The profiles where the LDS 
is found always have a shape similar to that in Fig. |3Kb) 
and Fig. |31 Intuitively, the existence of multiple locally- 
minimizing histories leading to the same p/ is due to 
the favorable action due to large a (p) on certain trajec- 
tories, utilizing densities on either side of the minimum 
in (J (p). A similar argument can be given for the ratio 
a (p) / D (p). The existence and exact location of the LDS 
depends on the full functional form of a (p) and D (p) . 
It would be of interest to find precise criteria. 



IV. STRUCTURE OF CUSP 

As discussed above, the structure of the LDS is similar 
in different models. Consider pf in some fixed 2d cross- 
section of the full phase-space, e.g., the cross-section de- 
fined in Eq. ([T]). As can be seen in Fig. EJb), the switch- 
ing line ends at a single profile (a point in the cross- 
section), which we denote by p™'*^(a;). This is much 
like a first-order transition line ending at a second order 
point. We now discuss the behavior of the LDF 
near p™*'^, as a function of p/ and N. As we now show, 
in the simplest scenario (j) behaves like in a Landau mean- 
field second order phase transition, or a "cusp catastro- 
phe" in the language of catastrophe theory [sO, IMj . 

The discussion builds on previous results pert aining to 
systems with few degrees of freedom [l^, ITtI - IIQI [33j . The 
singularity structure is well understood in such systems, 
where catastrophe theory is applicable. The extension to 
fields requires care, as we show below, see discussion at 
the end of this section. We first present the theoretical 
considerations. AppendiceslClandlDlverifv the prediction 
for the QS and BDI models. 




FIG. 6. Definition of quantities near the cusp. Dashed line 
switching line. 



at fixed p the minimum over J, subject to Eq. ([T}, is 
unique. It will therefore be convenient to work with the 
action after J has been minimized: 



s [p] = niin S [p, J] 



For a given pf on the switching line there are two his- 
tories, pi {x,t) and p2 {x,t), which minimize the action, 
as in Fig. [D We introduce 



Pf 



L/ r I iJbJu 



1/2 



as the distance of the final configuration pf from p™**^, 
see Fig. [6] We define a coordinate system, (a, 6) on the 
cross-section, with p™^^ at the origin, d directed along 
the switching line and positive on the switching line, and 
h orthogonal to the switching line. In analogy with Lan- 
dau mean-field theory, a plays the role of {Tc — T) and h 
the role of the magnetic field. 

Close to the cusp, when moving in the positive a direc- 
tion, for small enough s has two locally-minimizing 
solutions, pi and p2. Let (see Fig. ^ 

C')(^'^) = ^['°i(^'*) + ''2(a;,t)] , 
f 



and 



h(afi) {X, t) = 2 (x, t) - Pi {X, t)] , 
U(a,b) {X,t) = 5p(a,b)/ \\5p(aA)\\ ■ 



A = |l<5p|l 



(8) 



where we quantify the distance between two histories by 
\\5p\\^ = J [Sp (x, t)]^ dxdt. Here A plays the role of the 
amplitude of the order parameter. Note that at the cusp 
A = 0. 

On the switching line 6 = and s [pi] = s [p2] by defi- 
nition. Hence 



A. Multiple minima near the cusp 



Ha.b) (q) 



avq 



The action S [p, J] is a functional of both p and J. 
The dependence of S on the current J is quadratic, and 



has two minima, at Qmin = ±A. q is an "order param- 
eter" interpolating between pi and p2. Close to p™'*^ 
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the two minima should approach each other, and merge 
to a single minimum at p^^"^ . The simplest form which 
captures this behavior and is analytical in q, a and h is 



S(a,b) (q) = S(a,6) {q)~S(a,b) (0) c^q^ -ac2q'^ +Cibq , (9) 

with ci,C2,C4 > constants, s can also include higher 
powers of q, which would not affect the behavior at small 
a. At small a and b — 0, s (q) has two minima, at 



has at least one vanishing eigenvalue for the optimal path 
ending at p'f^. (This is a standard result in catastro- 
phe theory [3l|.) The corresponding eigenvector u {x,t) 
is precisely U(a,6) defined in Eq. ([8]) for p'j""^, i.e. 

u(x,t) = U(a^o+,b=o)- As s is minimal, H is positive 
semi-definite, and its entire spectrum is non-negative. 
We now assume that there is a single zero eigenvalue, 
followed by a finite gap. We then split the histories as 
follows 



± Va, hence A cx in direct analogy with Lan- p (j,^ ^ p^^^ (j,^ f- 5) + (^.^ ^) + p^ (3,^ f- p^ (^^ 5)) 



dau theory, with a mean-field exponent equal to 1/2 

In Appendix [D] we check the validity of Eq. ([9]) on 
the BDI model. We show that it indeed holds, but that 
higher order terms are still significant until close to the 



cusp point ( 



Pf 



cusp 

Pf 



10-2). For the QS model we 



use a different approach, see below. 



B. Effect of "soft modes" 

When Pf approaches p™^^ from the positive a direc- 
tion, the height of the action barrier separating the two 
local minima vanishes. This means that the contribution 
of the paths close to the minimal paths is enhanced. As 
we now show, this gives a new universal contribution to 
the probability of p^"''^, scaling logarithmically in N 



P 



cusp 

Pf 



exp —N(l) 



cusp 

pf 



+ \loeN 



(10) 



The universal factor 1/4 is known as the Berry index in 
catastrophe theory [30| . 

To see this, we go back to the path integral formula- 
tion, P[pf] — J DpDJexp{—NS[r,J]}. In this expres- 
sion, if the path integral is discretized then the measure 

is Dx = Yii (^VNdxi^ , where the y/N ensure normaliza- 
tion. For given p(x,t), S* is a quadratic functional in J, 
so J can be integrated out. For large N a saddle-point 
approximation gives 

In^y" DJexp{-NS[p,J]}^ ^ ~N mm S [p, J] = -Ns [p] 

and the path integral now reads 

P[pf] = J Dp{x,t)exp{~Ns[p]} . (11) 

Note that the correction to the saddle-point is N'^ in this 
case (3^ . 

We focus on final configurations in the cross-section. 
Since close to the cusp, when moving in the positive a 
direction, s [pf] has two solutions, this means that the 
Hessian matrix 



H = 



(12) 



Here pavg is defined as above when there are two minima, 
and is equal to the minimal history when it is unique, q 
is a numerical prefactor, and all other contributions are 
included in p^{x,t). Integrating out the p±_ directions 
we are left with an integral over q 

P[p/(a,6)]^e-^^(-^)(0)7Vi/2 f rf^exp {-iV5(a,fc) (g)} ■ 



where S(^a,b) (?) is defined in Eq. The N^^"^ comes 

from the definition of the path integral measure. The 
form in Eq. ^ was argued on the basis of the analyticity 
of s , justified by our assumption of the gap in H . The 
integral 

^ (TV, a, b) = TV^/^ ( dqexp [~N {cibq ~ C2aq^ + €49"*)] 



is known as the "cusp diffraction catastrophe" [s^l- We 
note two of its properties: (a) the "metastability" region, 
where the integrand has two local minima as a function 
of q is bounded by 6 cx ±0^^^. (b) tp {N,zi,Z2) has the 
scaling property (with u = N^/^q) 

ipN {a,b) = N^/'^ j dqexp [-N {cibq - C2aq^ + aq^)] 

^N^'^^(N^'^a,N^/%) , (13) 

where 5" (a, 13) — f dv exp [— (ci/Jw — 02^"^ + C4v'^)] has 
no N dependence. At (a, b) = (0, 0) this becomes 
i/;(0,0) = 7Vi/4^'(0,0). 

Therefore, at p'p''^ we have (j) p'J^''^ = S{a,b) (0)j ^^'^ 

P has an additional N^^'^ prefactor to the prob- 

ability distribution shown in Eq. (jlOp . This means that 
at the cusp the exponentiated N dependence has an ad- 
ditional non-analytic contribution, scaling as log N with 
a universal prefactor. The exponents 1/4, 1/2, 3/4 in Eq. 

were introduced in [34] and This implies that 
the logA^ correction in Eq. (|10p affects the probabil- 
ity in an elongated region of dimensions Aa x Ab ^ 
7V-1/2 X around p™"^ 

The above analysis relies on the assumption that the 
Hessian spectrum has a single zero mode followed by a 
gap. This can be generalized to situations where there 
is a finite number of zero modes followed by a gap, us- 
ing tools from catastrophe theory. In such cases, more 
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complicated singular structures will appear at p'p'^^, 
with modified universal exponents. The existence of a 
gap is expected to always hold in systems with finite- 
dimensional phase spaces. However, in the case of fields, 
where the phase space is infinite dimensional, the Hes- 
sian may be gapless. Then the analyticity of the action 
might fail altogether, as indeed happens in equilibrium 
critical phenomena [35l|. 

In Appendices ICl and [Pl we show that the assumptions 
indeed hold. Specifically, for the QS model, we prove 
that for specific types of profiles the assumption of ana- 
lyticity is justified. In addition, we calculate the Hessian 
spectrum numerically, and find a gap above a single zero 
mode. For the BDI model, we show numerically that the 
action indeed has a Landau mean- field form, Eq. 0. 




FIG. 7. (a) Toy model paths (thin lines) compared with cross- 
sections of exact paths for the SSEP model, (b) Same as (a) 
for BDI model, (c) Cusp area for toy vs. exact in the BDI 
model. BCs are pL = 0.2, pR = 0.9. 



V. THE CONNECTION WITH 
FINITE-DIMENSIONAL PHENOMENA - A TOY 
MODEL 

In an attempt to better understand the LDS in this 
system, we note that the minimizing histories p {x, t) in 
Fig. |4] appear to be quite smooth in x; this is sensible, as 
the field is constantly diffusing, making enduring, s harp 
spatial gradients improbable. This was studied in [13|. 
The smoothness motivates us to introduce toy models 
with a finite number of degrees of freedom, which cap- 
ture many of the essential features of the field models 
described above. 

To this end we discretize the field p {x,t), replac- 
ing it with a vector Pi{t), i = l,..^Np, correspond- 
ing to the density at the points Xi. Substituting Eq. 
([2]) into Eq. ([T]), the Langevin equation reads dtp — 

D (p) V p + \J a {p)r\ . For x G [0, 1], we take Xi — 

i/{Np -f 1), and obtain Np coupled Langevin equations 

dtPi = (Ax)"^ [A+1/2 {Pi+l - Pi) - A-l/2 {Pi - Pi-l)] 

+ (Aa;)"^ - V^i-i,i77i-i/2] , (14) 



where Aa:; = {Np + 1) ^, and po,pNp+i are assigned 
the boundary values /9b(0),pb(1) respectively, and 



-1, CTj.i+i are an ap- 



(»7i+i/2»yj+i/2) = (Ax) iV i(5ij. A 
propriate choice for D (p) , a (p) for x between Xi and 
Xi+i. We choose A,i+i ^ ^[D {pi) + D (pi+i)]. A sim- 
ilar choice can be made for (7^ ^^i, but we use (Xj i+i = 

2£)i,,+i (p,+i - p,) [/' (pi+i) - f ip^)]~\ where /' (p) is 
given in Eq. This has the advantage that when 

the BCs are equal, po — PNp+i, the system of Langevin 
equations satisfies detailed-balance |38| with respect to 
the potential (p{{pi}) = AxJ2if{Pi)- The discrete 
Langevin equation converges at high Np to the field- 
theory (as can be seen by writing the action) . As the his- 
tories are smooth we expect rapid convergence for long 
wave lengths. We therefore use the lowest non-trivial 



discretization, Nr, 



2, which can accommodate non- 



equilibrium phenomena. In this case pi,2 correspond to 
"coarse-grained" densities at x = 1/3,2/3 respectively. 

The minimizing histories for the toy model with Np = 
2 can be obtained using standard techniques from low- 
noise finite-dimensional systems. Using the approach 
outlined in the introduction, a set of coupled ordinary 
differential equation is obtained, whose solutions are the 
extremizing histories. The equations are solved numeri- 
cally using a shooting method [l4|. 

In Fig. [7] we check the above approach on the sim- 
ple symmetric exclusion model (SSEP) with D — I and 
a = 2p(l — p), which does not feature a cusp, and for 
the BDI model. We plot trajectories of the toy ver- 
sion of the SSEP (Fig. [T^a)) and the BDI model (Fig. 
WO^)) in the (pi,p2) plane, against the {pi/s, P2/3) tra- 
jectories of the full models. In addition, the metasta- 
bility region for the toy model and in the cross-section 
of the exact dynamics are plotted. The qualitative pic- 
ture is similar - a metastability region appears in the 
quadrant pi/3 > 0,p2/3 < 0, at approximately the 
same location as in the exact field solution. This is 
expected to have a close relation to the breaking of 
detailed-balance, which is easy to visualize in the two- 
dimensional toy model. Define the two- variable Langevin 
equation dxi/dt = A'i(xi,X2) -|- 2 ^y'^i' ^^'^ "^^^ 

Q = BB"^. Detailed-balance is satisfied if V0 — Q^^K, 
or V X (Q^^K) = 0. Therefore, in two (phase-space) 
dimensions w = V x (Q~^K) is a scalar which quan- 
tifies the breaking of detailed balance. It is shown in 
Fig. [5] Minimizing trajectories passing through regions 
with a; > (w < 0) bend counter-clockwise (clockwise). 
Therefore, a gradient in u can cause trajectories to focus 
and cross, creating a cusp. A similar picture has been 
discussed for other finite-dimensional systems jTsf . 

In summary, the above analysis shows that much of 
the phenomena observed can be captured by simplified 
finite-dimensional models, making concrete connections 
to previous works on such models. 
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CO 



1/3 



FIG. 8. Map of the measure for breaking of detailed- 
balance, in the BDI model. Clockwise currents for positive 
uj. In black: selected baths. 



VI. DISCUSSION 

Many interesting questions remain to be studied. First, 
it would be of interest to find precise conditions or bounds 
for the occurrence of the LDSs, depending on the bound- 
ary conditions and model parameters. In addition, a sim- 
ple picture for the mechanism leading to the existence of 
multiple histories ending at the same final profile is lack- 
ing. The mechanism suggested by the current authors in 
[2]| was flawed [36j . 

The singularity described above has the simplest possi- 
ble structure. More complicated singularities are in prin- 
ciple possible. For instance, in catastrophe theory, richer 
structures have been analyzed. Observing them in full 
requires one to look at higher dimensional cross-sections. 
It would be interesting to find which of them exist in dif- 
fusive models, and for which models. Even more exciting 
possibilities exist: Fields allow for the possibility that 
the Hessian discussed in Sec. llVl is gapless. Thus, it may 
even be possible to find singularities that are beyond the 
realm of catastrophe theory. 

Finally, it would be very interesting to look at these 
LDSs in higher dimensional systems. This is now possible 
using the numerical technique described in [T3| , and used 
in the present paper. 
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Appendix A: Calculating D (p) and a (p) for the 
driven Ising model 

As shown in [2I, l37l|, for each parameter set {e,6) one 
can write implicit analytic equations for D (p) which 
can then be inverted numerically. Then a (p) is ob- 
tained via the fiuctuation-dissipation relation, a (p) = 



2kBTp'^K (p) D (p) where k{p) is the compressibility 
For equilibrium BCs this model admits an Ising measure. 
To find D (p), we use the relation [s^ 



X 

1^) (related to the compress- 



where x = E< (("•i"o) - P^) 
ibility k (p) by x = ksTp^n (p)), ji^i+i is the current 
(number of particles per unit time) from site i to site 
i + The averages are taken with respect to the equilib- 
rium probability distribution, and {ji^i+i) = due 
to the symmetries in equilibrium. One then finds a using 

(7 = 2xi? = 2 0,,+i) . 

To calculate (ji^i+i) note that as the rates depend on 
the four sites around a bond, we have that 

{J^,^+l) = (1 + <5) PolOO + (1 + ^1100 
+ (1 - £) Poioi + (1 - <5) PllOl 

where Pqioo is the probability of the pattern 0100, and 
similarly for others. Using the transfer-matrix technique 
[S^], one can calculate these probabilities and obtain 



\[\ + 5{l-2p)]-e^Ap{l- p) 
A3 



where 
A = 



V4p(l-p) V4p(l-p) 



l + £ 



1/2 



It remains to find p and k, which are both given in terms 
of (3 (recah that exp (4/3) = (1 -t- e) / (1 - e)), and h = 
/3p: 



sinh h 




Ve4^+sini?7i, 
e^^ cosh h 



4 (e4/5 -I- sinh^ h) 



3/2 



In order to obtain D{p),a{p), we calculate 
a {h) , D {h) and p [h) for a wide range of /i, and nu- 
merically invert the last to find D [p) — D{h{p)) and 
CT (p) — a{h{p)). Fig. [2](a) shows D (p) and cr (p) for 
{e,S) = (0.05,0.995). 

As a check, we note that for the simple symmetric ex- 
clusion process Ij 5 ^ e = 0, and one finds 



p = - (1 + tanh h) 



1 



4 cosh h 



and 

A = 2 cosh /i ; {jt.i+i) = A"^ 
so that D — 1 and cr = 2p (1 — p) [Ij. 



1 



4 cosh h 
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P(x) 

■ - -g,M 

■ - -g2('<) 

■ - -g,(x) 



FIG. 9. Density profile p{x) of tfie step form (Eq. (IBip ). and 
three g{x) solutions. Here pA — 4:, ps = — 5,p_i = — 2,pi — 
3. 



Appendix B: Existence of multiple extremal 
solutions in the QS model 

In this Appendix we prove that for the QS model, 
which has D — 1 and <J {p) — + 1, and for any non- 
equilibrium BCs, there exists a LDS for some profiles. 
Here it will be far more convenient to work in the do- 
main X G [—1,1]. The results in the new domain are 
simply related to the results in the original domain (40| . 

The BCs to Eq. ^ are denoted by p_i = pL and 
P+i = PR- 
Claim 1 For any BCs p_i =/= p+i, there exists a profile 
Pf {x) for which Eq. (0j is satisfied by more than one 
solution with g (±1) — p±i. 

Proof. Using the symmetries p —p and x ~x it is 
enough to consider the case p-i < pi, and < pi. 

We proceed by an explicit construction of p/. That is, 
given p±i we construct a function pf {x) for which Eq. 
([6]) is satisfied by more than one function g [x) , which 
also satisfies the boundary-conditions. The profile pf (x) 
will be a piecewise-constant function composed of two 
flat regions, of the form 



Pf (x) 



PA -I < X <0 
Pb < X < 1 



(Bl) 



where pa,Pb are (constant) numbers which specify 
Pf {x), see Fig. [S] Note that pj {x) does not have to 
be continuous, nor to satisfy the BCs, hence pa, Pb are 
not restricted in any way. The solutions g (x) will be put 
together by solving Eq. ([6]) for x < and x > sepa- 
rately (each with its corresponding boundary condition) , 
and matching the solutions by demanding that g (x) and 
g' (x) are continuous at x = 0. 

For a region with constant p (x) = p, and given g (xi), 
Eq. has an (implicit) analytic solution 



g{x) gpatani> 



:dlp — Cl {x — Xl) 



(B2) 



where Ci is a free constant. Differentiating both sides 
with respect to x we find 



atan g{x) 



Cl 



(B3) 



'l + gixY 
and using Eq. (jB2p for ci, g' (x) reads 

\/^ + 9{xf 1 



/(x) 



X-Xl Jg(^,) ^/l+W' 



:dlP . (B4) 



Note that ci no longer appears in this equation. Instead, 
this is a relation between g (x) and g' (x) . Let gA (x) be 
the solution given in Eq. (jB2[) with xi = —1, g (xi) — 
p_i and p = pa: 



9Aix) atanV 



d^/j = CA (x + 1) 



(B5) 



for — 1 < X < 0. This defines a one-parameter family of 
solutions, according to the value of ca- Similarly, gs (x) 
is defined by 



rgi 
J PI 



Sb(^) gPB atani/> 



dTp = cb {x — 1) 



(B6) 



for < X < 1. Any solution of Eq. © with p (x) of the 
step form defined in Eq. (jBip is composed of solutions 
9A (x) , gB (x) satisfying gA (0) = 5s (0) and g'^ (0) = 
g'g (0). The derivatives at x = are given by 



9'a (0) 



9'b (0) 



I + gA (0) r9A{0) atan V 



gPA atan gA(0) 



1+5B (0)^ 



gP_B atangA(O) 



p_i yjl + i'- 
Pi gPs atan ip 
b(0) \/1 + V'^ 



di) 



(B7) 



We note that: 



(a) g' (x) does not change sign. As we are interested in 
solutions with p_i < pi, we only need to consider 
solutions with g' (x) > 0. 

(b) From (a) it follows that p_i < gA (0) = gB (0) < 
Pi- 

(c) It also follows that if gA (0) = p-i then (0) = 0, 
and if gB (0) = pi then g'^ (0) = 0. Similarly, if 
9A (0) = pi then g\ (0) > 0, and if gB (0) = P-i 
then g's (0) > 0. 

Consider now (0) and g'^ (0) as a function of g (0). 
A solution g (x) on the entire segment [—1, 1] is obtained 
when g'^ (0) = g'^ (0) for the same g (0). Remark (c) en- 
sures that they cross at least once; But they may also 
cross more than once, see Fig. [TUl The number of cross- 
ings depends on pa,Pb- We will show that there always 
exist patPb for which the graphs cross more than once. 
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FIG. 10. g'j^ (0) and g'g (0) plotted as functions of g (0) for two 
functions p{x). In the upper pane the graphs cross only once, 
indicating a single g (x)-solution. In the lower pane, done 
with parameters of Fig. (|9]), they cross three times, resulting 
in three different p-solutions. (Note that the g-values at these 
crossings indeed correspond to g (0) of the solutions in Fig. 
Q). BCs for both panels are p_i — — 3, p+i = 5. Upper 
panel: pA — S, pB = —2, lower panel: pA — 4, ps = —5. 
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FIG. 11. g'A [g (0) ; pA = 100; p-i] for different p_i. Gray line: 
large p-expansion, Eq. 



Motivated by the fact that the cusp singularities always 
appear at the lower right corner of our phase-space cross- 
sections, we consider the limit where pA is a large positive 
number, and ps — —pA- 

Denote by 5^ [5 (0) ; p^; P-i] the value of 5^(0) 
as a function of g (0) , pA and p-i , and similarly 
g'B [9 (0) ; pb', Pi]- The analysis which follows is done 
for g'j^ (0); similar results are obtained for g'^ (0) since 
9'b [5(0) ;pb;Pi] = g'A [-5(0);-Pb;-Pi]- To better un- 
derstand 5^ [g (0) ; Pa; P-i] at large pA, we plot g'^ [g (0)] 
for PA — 100 and different p_i values, see Fig. [TTJ As can 
be seen, the different graphs rise quickly from g'j^ (0) = 
at g (0) = P-i, and join a common function. This is 
formulated by the following Lemma: 

Lemma 2 Expanding around pA ^ 00, we have for any 
g(0) >p_i 

g'^[g-p^;p^,]^l±^^ii^±^ + 0{p-/) . (B8) 



Proof. Rewrite Eq. dBT]) for g'^ (0) as 

p-^PA [atan g—atan ip] 



9'a (0) = 



where here and in the next equation g stands for gA (0). 
For p^ — >■ 00 a saddle-point approximation can be pre- 
formed. The exponent —p^i [atan — atan -0] is domi- 
nated by small values of g — ip, i.e. close to the upper 
bound of the integral, and atan<^ — atan-i/; can be ex- 
panded to second order 



atan g — atan ^ = 



g-i' , ig-->P) 



1 



ii + g'Y 



+ (g-^r 



r /2 

In addition, the denominator (l + ■(/;^) is expanded 
to second order around g. The resulting expression in- 
volves Gaussian integrals which can be integrated, with 
the lower integration limit set to — cx). Finally, we expand 
the result (containing error-functions, etc.) to second or- 
der in 1/pa around p^i — >■ 00, and obtain Eq. (IB8|) . ■ 

The expression in Eq. (jB8|) does not depend on p_i, 
as expected from the reasoning alluding to Fig. [TT] Sim- 
ilarly, for g'g [g (0) ; ps; pi] we have, for pB —00, 



g's [g;pB;pi] = - 



1 



5(1 



PB 



Pi 



o 



[pb') 



(B9) 



We are now in a position to construct p/ (x) with three 
solution to g: given the BCs, choose some crossing value 
gc G (p-iiPi) (this is where the non-equilibrium con- 
dition p_i 7^ pi enters). For a given pA the condi- 
tion g'j^ [gc, PA] P-i] = g's [gc, Pb\ Pi], together with Eqs. 
(|B9|) reads, to second order in p^^, p^^. 



1 



gc 



gc 



{i + gl) i + gl , gc{i + gl) 



PA 



P\ 



PB 



Pi 



or 



^fi-^U-^ 

PA V PA J PB 

Solving this for pB find 



1 - 



gc_ 

PB 



(BIO) 



(Bll) 



-1 



PB 



1 



PA 



PA 



2 1- 



(B12) 



PA 



Pa 



Fixing gc, the expansions Eqs. (|B8|) . (|B9|) guarantee that 
as PA grows, with pB given by Eq. (|B12p . a crossing point 
will appear in a neighborhood of gc, and approach gc for 
Pa ^ 00. The root of the quadratic equation leading 
to Eq. (jB12[) was chosen so that g'A,g'B '^i^l cross from 
9'a [9] > 9'b {9] for 9 < gc, to g'^ [g] < g'g [g] for g > 
gc, see Fig. [TH This, together with g'^ [p_i;pA;p-i] = 
and g'g [pi;pB;pi] — 0, guarantees that the functions 
5^ [5;Pa;P-i] and 5b [g; ps; pi], wiU have two crossing 
points in addition to the crossing point near g = gc- For 
large pA these will be at values of g close to p_i and p+i, 
see Fig. ^ 
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FIG. 12. Constructing a solution. Choose a crossing value 
Qc (here (jc = —1) and some large pA- pB is fixed so that 
the crossing is approximately at Qc (Upper pane). The full 
solutions will feature this crossing with two more crossings, 
close to the boundaries. 



Appendix C: Cusp Structure and Hessian spectrum 
in the QS model 

In this appendix we study the cusp structure, and the 
spectrum of the Hessian matrix H at p™*^, for the QS 
model. First, we prove that for the QS model with pro- 
files defined as in Appendix [B] the framework of catas- 
trophe theory is applicable. More precisely, there exists 
an analytical function F (pa, Pb) on the two-dimensional 
cross-section parametrized by {pA, Pb) as in Appendix iBl 
and for which every extremum of F corresponds to a sin- 
gle extremal history leading to Pf (pa, Pb), as defined in 
Eq. (|BT]) . 

To construct the function F, we note that Eqs. (jB7[) 
give us analytical expressions 3^(0) = fAigo,PA,PB) 
and g'siO) = /b (go, PA, Pb)- We drop the p-i,pi de- 
pendence, which are kept fixed. Let 

F {go; PA, pb) = Wa (0) - g'B iO)f ■ 

F {(j)Q; PA, Pb) is analytic in all its variables, and 
F (go; PA, Pb) — iff the solution g is an extremal solu- 
tion. Moreover, in the vicinity of the cusp, dFjdg^ = 
iff F (f/o; Pa,Pb) = 0, i.e. is a local minimum as a func- 
tion of go- Therefore F acts as a "gradient map" (in 
the sense of Catastrophe Theory), with pa,Pb the con- 
trol variables, and go the state variable. Accordingly, the 
cusp structure (regions in {pa,Pb) where F (go) = has 
two solutions) is expected to be mean-field. 

To show how F is used, we briefly review the argu- 
ment for the cusp structure, which is essentially a Lan- 
dau mean-field argument. The cusp point is a special 
point {p'^A^^, Ps"*^) where at the minimal go, d'^F/dg^ — 
= d^F/dgQ hold (the two conditions explain why it is a 
point, or a set of isolated points, in the (pA, Pb) plane). 



Let xa,b = Pa,b — p'a"^ ■ In the vicinity of the cusp, 
expand F to forth-order [here the analyticity is crucial), 
F — agQ + l3gQ + jg^ + S, where a,/3,7, 5 depend on 
XA,B- We assume that a ^ 0; vanishing a would be non- 
generic, i.e., could be remedied by a small change in any 
additional parameters, such as the BCs p^i,pi or the 
noise function a. Then a local change of variables can 
be performed to bring F to the form F = ^g^ -\- ag^ -\- b, 
where at the cusp a ~ — b. The region where F has 
two local minima is bounded by 6 = i^^a^/^. Near 
the cusp a,b can be expanded to first order in xa,b, so 
the power-law relation between b cx ±a^/^ will apply to 
a rotated frame of xa b ■ 



1. Spectrum 

As discussed in the main text, for the "Landau mean- 
field" , catastrophe theory to hold, one must have a gap 
in the spectrum of the Hessian H. This can be tested 
numerically, by evaluating the action S for the extremal 
solution p (x, t) leading to pf {x), and calculating the Hes- 
sian, Eq. p^ . by varying p jointly at pairs of points 
(xi,ti) and {x2,t2). The eigenvalues of H can then be 
calculated for different pf profiles close to p/ = Pcusp- 

To calculate H one thus needs to locate Pcusp- This is 
most easily done for the QS model with BCs p^ = —pR, 
where Pcusp of the form of Eq. ([7]) must have ai = 0. Fig. 
[T^ shows the bottom of the spectra of H for different p/, 
starting from p and ending at Pcusp- One can clearly 
see a single eigenvalue going to zero, in agreement with 
the analysis in the paper. The rest of the eigenvalues 
remain away from zero, without closing the gap. This 
validates the analysis carried out in the main text for the 
QS model. 



Appendix D: Cusp structure in the BDI model 

In this Appendix we check the validity of Eq. ([9|), 
which predicts the structure of the cusp. In the BDI 
model the diagonalization of the Hessian H gave incon- 
clusive results. We suspect that this is due to the diffi- 
culty of locating pcusp with high precision in this model. 
As we show now, the predictions of Sec. IIVI hold only 



cusp 

Pf 



< 10- 



very close to p™"'', when pf 

To compare Eq. ^ with numerics, it is more con- 
venient to use a different form, which does not re- 
quire knowing the precise position of p™"^. Noting that 
S(a,o) (A) = s [pi] and S(^a,o) = « [P2], we find that 

A oc y/a. Therefore we expect that for 6 = 



(Dl) 



(A) - s,=o (A) « ( -y 



where y — <z/A. 
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FIG. 15. (a)ThefunctionS'g(A)-5'g=o(A)/A'' as a function 
of 2/ at a = 0.012 (solid line). A clear deviation is seen from 
a fit to a quartic function (dashed line), (b) The function 
Sq (A) — 5*9=0 (A) / A* for different values of a (dashed lines). 
Fitting the functions to a polynomial of power 8, and plotting 
only the quartic part, the collapse improves significantly (solid 
lines) . 



FIG. 13. Bottom of Hessian spectrum for the QS model at dif- 
ferent p/ profiles, equaly spaced between p (leftmost) to p^""'' 
(rightmost). BCs are ph ~ —3,pR = 3. A single eigenvalue 
approaches zero, while the gap above it is maintained. 
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FIG. 16. Fits of C4. (a) C4 as a function of A^. (a) a as 
a function of a. Guidelines (dashed lines) represent the ex- 
pected slope of the functions. Extracting a from fits that 
also include higher powers one obtains data that better fits 
the expected slope. 



FIG. 14. Pairs of locally minimizing histories, leading to 
points on the switching line (solid and dashed lines). The 
history leading to p^"^ (circle) is also plotted (bold line). 



As an example we consider the boundary-driven Ising 
model, with {e, 5) = (0.05, 0.995) and pL = 0.2, pn = 0.8. 
Examples of pairs of locally minimizing histories leading 
to configurations on the switching line are shown in Fig. 

M 

Fig. ITSl a) shows the function [sq (A) — Sq=o (A)] /A*^ 
as a function ofy at a ~ 10"^, together with a quartic fit, 
which shows clear deviations from this form. This means 
that even at this distance Pf — p^"^ 10~^ there are 
significant contributions of higher powers to Sq. Due to 
these higher powers plotting [sq (A) — Sq^o (A)] / vs. 
y for difi^erent a values in the range 5 • 10~^ < a < 0.2 
does not collapse the data as expected. We therefore fit 
the functions Sq (A) — Sq=o (A) to polynomials of order 
four and higher, and plot C4, the prefactor of q"^, as a 



function of A^, see Fig. ITSf b). The expected power- 
law is not obtained for a quartic fit, but improves when 
the fits include higher order terms, see Fig. [TCTa). This 
means that C4 in Eq. ^ can indeed be taken to be 
constant. Finally, one can also fit Sq (A) — Sq^o (A) oc 
'^^ {\y^ + fy^)' by fitting the position of p™"^, see Fig. 

The presence of strong higher powers as close as a ~ 
10~^, see e.g. Fig. IWb). can be understood as follows: 
the non-linear terms come from the different action at 
the two paths pi and p2- As the distance A between the, 
scales as A (X ^/a, pi— P2 at some space time point, can be 
of order 10~^ even for a ^ 10^^, and the two minimizing 
paths can see very different behavior of D (p) ,a (p) along 
the paths. This sensitivity explains why computing H 
directly is difficult: one needs A to be small, but as o oc 
A^ , one needs the distance a from p™*^ to be very small. 

From the above we conclude that at the cusp there 
is a soft mode, a direction along which the minimiz- 
ing history p™"? (x, t) has a zero second derivative: 
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(Ps [p^'^^P {x,t) + qUa=o {x,t)] /dq^ ~ [(Psq/dq'^'^^_^ = 0. This is indirect evidence that H has at least one vanish- 
ing eigenvalue. 
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